Green's Function Approach to the Bose-Hubbard Model 
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We use a diagrammatic hopping expansion to calculate finite-temperature Green functions of 
the Bose-Hubbard model which describes bosons in an optical lattice. This technique allows for a 
summation of subsets of diagrams, so the divergence of the Green function leads to non-perturbative 
results for the boundary between the superfluid and the Mott phase for finite temperatures. Whereas 
the first-order calculation reproduces the seminal mean-field result, the second order goes beyond 
and shifts the phase boundary in the immediate vicinity of the critical parameters determined by 
the latest high- precision Monte-Carlo simulations of the Bose-Hubbard model. In addition, our 
Green's function approach allows for calculating the excitation spectrum at finite temperature and 
for determining the effective masses of particles and holes. 



Introduction: Ultracold bosonic gases trapped in the 
periodic potential of optical lattices represent tunable 
model systems for studying the physics of quantum 
phase transitions 0, H, H, |4j. They are described by 
the Bose-Hubbard Hamiltonian which decomposes into 
two parts: H = Ho + H\. The first term Hq = 
Si [Ufii(ni — l)/2 — [ifii], with the on-site energy U and 
the chemical potential fi, describes the repulsion of more 
than one boson residing on a lattice site. It is local 
and diagonalizable in the occupation number basis for 
any lattice site. The second term Hi — —Jij : a\dj 
with the hopping matrix element Ji j — J if the lattice 
sites i and j are nearest neighbors and J — other- 
wise describes the hopping between two sites due to the 
quantum- mechanical tunneling effect. The competition 
between the two energy scales U and J determines the 
existence of two different phases. When the on-site en- 
ergy is small compared to the hopping amplitude, the 
ground state is superfluid (SF) as the bosons are delo- 
calized in a phase coherent way over the whole lattice. 
In the opposite case, where the on-site interaction domi- 
nates over the hopping term, the ground state is a Mott 
insulator (MI) where each boson is trapped in one of the 
respective potential minima. 

This characteristic quantum phase transition of the Bose- 
Hubbard model has been studied extensively both with 
analytical |, S 0, H H and numerical [H El El 
methods for zero temperature, while less literature ex- 
ists on the finite-temperature properties of this transition 
0, EH EH- In this letter we work out a Green's func- 
tion approach in order to determine the MI-SF phase 
boundary and the excitation spectrum in the Mott phase 
beyond the mean-field level for finite temperature with a 
hopping expansion. Our analytical findings compare well 
with the latest high-precision Quantum Monte Carlo sim- 
ulations and allow to propose a thermometer for bosons 
in optical lattices. In the following we restrict ourselves 
to a spatially homogeneous system and neglect the effects 
arising from the additional harmonic confining potential, 
which is present in all experimental settings, like the for- 



mation of a shell structure [l6(. However, these effects 
could be taken into account by applying the local density 
approximation where the external potential is taken into 
account in form of a spatially dependent chemical poten- 
tial. 

Green's Function Approach: As all single-particle proper- 
ties of a quantum many-body system are contained in its 
Green function, we base our calculation on this quantity. 
Because we are interested in describing a system at non- 
zero temperature, we use the imaginary-time formalism 
[13, Ell- Therein, the single-particle Green function is de- 
fined as the thermal average of the time-ordered product 
of a creation and an annihilation operator in Heisenberg 
representation 
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with H = 1, /3 = 1/fcsT, and the partition function 
Z = Tr{e~@ H }. Because it is not possible to obtain ana- 
lytic expressions for the eigenstates and eigenenergies of 
the full Bose-Hubbard Hamiltonian, we cannot calculate 
the Green function exactly. Instead, we aim at a pertur- 
bative treatment and calculate this quantity as a power 
series in the hopping matrix element Jij. As that pa- 
rameter is small for the Mott phase, where the lattices 
are deep and the interaction between particles is strong, 
we refer to this treatment as a strong-coupling expan- 
sion. In order to employ this perturbative expansion for 
finite temperature, we make use of the Dirac interaction 
picture and write the imaginary-time evolution operator 
in form of a Dyson series, 



U(t,t q ) = Texp 
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where the time dependence of the Dirac-picture operators 
is determined by the local Hamiltonian H - With its help 
we can write the Green function as 
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where the time-ordering operator acts also on the 
time variables which are resulting from the expan- 
sion of the Dirac imaginary-time evolution operator in 
(|2"|). When we now expand perturbatively in pow- 
ers of the tunneling matrix element J, the (n — l)th 

order contribution G$" ^ in ([3]) turns out to de- 
pend on the n-particle Green function of the unper- 



turbed system as Gn^ (t- 
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In order to de- 



termine Gn we cannot use standard Wick's theorem be- 
cause the unperturbed Hamiltonian Hq is not quadratic 
in the Bose operators. Although the lack of Wick's the- 
orem in the present situation prevents us from using the 
powerful perturbative technique based on standard Feyn- 
man diagrams, we can nevertheless simplify G^ by de- 
composing it into cumulants. To this end, we follow 
an approach reviewed by Metzner 19] in the context 
of the fermionic Hubbard model which describes elec- 
trons in a conductor. This decomposition is based on 
the important observation that the on-site Hamiltonian 
Hq is local. Consequently, the unperturbed Green func- 
tions G^ are also local and can be decomposed into 
time-dependent cumulants . For instance, we have 
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where = ^ 



-0E n 



is the unperturbed partition 



function for a single-site system with the on-site energy 
eigenvalues E n = Un(n — l)/2 — fin. With this decom- 

(n) 

position, we can represent G\ diagrammatically: We 
denote an n-particle cumulant at a lattice site by a ver- 
tex with n entering and n leaving lines with imaginary- 
time variables, so we have, for instance, for the first two 
cumulants 



C[°\r'\r) 
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Furthermore, the hopping matrix element is symbolized 
by a line connecting two vertices: 



Ji- 



(7) 



With all this, we can set up the diagrammatic rules for 
calculating the nth order contribution of the Green func- 
tion in J. First: Draw all possible combinations of ver- 
tices with total n internal as well as one entering and one 
leaving line. Second: Connect them in all possibles ways 
and assign time variables and hopping matrix elements 
to the lines. Third: Sum all site indices over all internal 



lattice sites and integrate all internal time variables from 
to (5. 

We also note here that we have to sum all site indices over 
the whole lattice, no matter whether two sites in a dia- 
gram coincide or not. We make use of the translational 
invariance in imaginary time and transform all expres- 
sions to Matsubara space. In the second diagrammatic 
rule the integrals over the time variables have then to be 
replaced by sums over all bosonic Matsubara frequencies 
uj m = 2um/ [3 with integer m where the sum of the in- 
coming frequencies must equal the sum of the outgoing 
ones. 

The formalism developed so far allows for calculating the 
Green function to any given order in the tunneling ma- 
trix element J. But because one of our main goals is to 
describe the phase transition between the Mott insulator 
and the superfluid phase, and it is well known in the the- 
ory of critical phenomena that such a transition is char- 
acterized by diverging long-range correlations [rsl . [20| . 
we must employ a non-perturbative method which is 
achieved by summing an infinite subset of diagrams. In 
order to perform this task, we introduce the sum of all 
one-particle irreducible diagrams including their respec- 
tive symmetry factors and multiplicities as 



Ci(w m ,k) = 
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The self-energy which describes the movement of a single 
particle in a many-body environment is defined as S = 
1/Gi —1/Gi [17]. For our case, it can be written in the 
formS(^ m ,k) = l/Gl 0) (w m )-l/CiK l ,k) + J(k), where 
J(k) = 2J^^ =1 cos(fc„a) is the Z?-dimensional lattice 
dispersion. For instance, the first order in the self-energy, 
which corresponds to a summation of all simple chain 

diagrams, yields Gi(w m ,k) = l/G^ (co m ) — J(k) 

with Gf\uj m ) = J2 n [( n + l )/( E n+i - E n + iuj m )- 
n/(E n - + iuj m )]e~P E "-/Z(°\ We note that per- 

forming the limit r' \ r in the Green function allows 
to obtain the quasi-momentum distribution needed to 
explain experimental time-of-fiight absorption pictures 

Hi mil. 

Phase Boundary: The boundary of the phase transi- 
tion is characterized by diverging long-range correlations 
[HI, H(| ■ Thus, we set k = and solve for the value of J, 
where the Green function diverges, which is only possible 
for u> m = 0. This yields up to first order in J 

2DJ r = M <- , (9) 



n+l 



which coincides with the finite-temperature mean-field 
result 0, EE HI- The zero-temperature limit of ©, 

i.e. 2DJ C = 1/ ( E Jl+l En - g „_g„_i )' agrees with the 
seminal mean- field result of Ref. [l[. The reason for 
this agreement is that each approximation becomes ex- 
act in the limit of infinite spatial dimension. In order to 
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FIG. 1: Quantum phase diagram for D = 3. Thick (solid): 
Second (first) order, JcbT = 0.1(7. Dashed (dotted): Second 
(first) order, T = 0. Dots: QMC Data for T = [ill. 



FIG. 2: Dispersion relation of the pair excitations for T = 0. 
Thick lines: Second order. Thin lines: First order. Solid: 
J = 0.029 (/. Dashed: J = 0.025 17. 



see that one must suitably scale the hopping parameter 
[H US H3- When we define J = 2DJ, the contribu- 
tion of the fcth order chain diagram is proportional to J k 
because there exist 2D possibilities in a chain diagram 
for every hopping line to connect to neighboring sites. 
The lowest-order term neglected by that summation is 
the one-loop diagram in (jSJl which has two internal lines 
but only one free index and is, therefore, proportional 
to 2DJ 2 = J 2 /(2D). Thus, it vanishes in the limit of 
D -> oo. 

Taking now the one-loop diagram into account yields 



1 
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where S( 2 )(o; m ,k) = Gf B) (w m )/ ]Gf\u m ) 
second-order self-energy with the value of the one-loop 
diagram G\ (ui m ). The analytic formula for the phase 
boundary resulting from this formula is shown in Fig. [TJ 
This result coincides in the limit of T — > with the effec- 
tive potential formalism employed in Ref. [9( ■ We empha- 
size that, unlike the first order the one- loop corrected 
result depends on the system dimension in a non-trivial 
way. Furthermore, the corrections are larger in two than 
in three dimension which is consistent with the fact that 
the first-order results become exact in the limit D — > oo. 
For finite temperature, the phase boundary is shifted to- 
wards larger values of J c as thermal fluctuations suppress 
quantum correlations which are responsible for the for- 
mation of the superfluid. This effect, which occurs both 
in first and in second order, is most important between 
the Mott lobes as fluctuations are strongest when the av- 
erage particle number is not near an integer value. The 
correction to the mean-field result arising from the one- 
loop diagram, visible in the difference between first and 
second order curve, plays an important role only near the 
tip of the Mott lobe. This feature stems from the fact 
that quantum fluctuation are notably increased when the 
system approaches the quantum critical point at the tip 
of the lobe [4|. Thus, we can say the quantum phase dia- 
gram consists of a thermally dominated region where the 



influence of thermal fluctuations is large and a quantum 
dominated region where the quantum corrections from 
the one-loop diagram are most important. 

Excitation Spectrum: For a system in the Mott phase 
three different types of excitations exist. 1. The addition 
of a particle from the environment (particle excitation). 
2. The removal of a particle (hole excitation). 3. The 
creation of a particle- hole pair (pair excitation) . The last 
one is important from a physical point of view because 
it is also possible in an isolated system as realized in 
the current experiments. In the superfluid phase, there 
are additional excitations corresponding to fluctuations 
of the phase of the superfluid order parameter. However, 
these features can not be investigated within the present 
formalism but with the related effective action method 

In order to obtain the spectrum of the quasi-particles, 
which is given by the poles of the real-time Green func- 
tion for finite temperature we must analytically continue 
our imaginary-time result. This is achieved by perform- 
ing the replacement iui m — > ui + ir/ with rj — > [17] . The 
first and second order, respectively, yield excitation spec- 
tra of which the former one agrees in the limit T — > with 
the mean-field result from Ref. Q ■ In the following we re- 
strict ourselves to the most interesting case D — 3. Both 
finite temperature and one-loop corrections are most ef- 
fective for small wave numbers, the former because the 
effect of temperature is the suppression of quantum cor- 
relations which mainly exist for long wavelengths, the 
latter because the dominant fluctuations near a quantum 
critical point are the ones with vanishing wave number 
as shown in Fig. [5] 

All spectra show a characteristic energy gap which van- 
ishes at the critical point, i.e. at the value of J c at the 
tip of a Mott lobe. Our prediction for the temperature 
dependence of this energy gap is shown in Fig. [3] As it is 
an experimentally accessible quantity [ll| it could serve 
for determining the temperature of bosons in an optical 
lattice. Furthermore, we determine for both the quasi- 
particles and the quasi-holes an effective mass which is 
shown in Fig. [U Note that they become massless at 
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FIG. 3: Temperature dependence of gap for unit filling. Solid: 
J = 0.008 U. Dashed: J = 0.012 U. Dotted: J = 0.025 U. 




FIG. 4: Effective masses of quasi-particles (solid lines) and 
quasi-holes (dashed lines) for T = 0. Thick lines: Second 
order. Thin lines: First order. Dots: QMC data [ill ]. 



the critical point as a consequence of the U(l) symmetry 
breaking at the second-order phase transition. 

Conclusion and Outlook: We have presented a pow- 
erful formalism to calculate the Green function for the 
Bose-Hubbard model in the Mott phase. It allowed us 
to improve the mean-field phase boundary in an analytic 
way for finite temperature. At absolute zero our result 
deviates from the recent Quantum Monte-Carlo studies 
of Ref. by only 3% for D = 3. For finite tempera- 
ture first analytic results beyond mean-field theory have 
been presented and the importance of both thermal and 
quantum fluctuations in the different regions of the phase 
diagram has been clarified. In addition, we have calcu- 
lated the excitation spectrum of the quasi-particles and 
derived its effective masses and compared them to numer- 
ical findings. We have also investigated the characteristic 
energy gap and determined its temperature dependence. 
In order to obtain more accurate results, it would be 
desirable to extend the hopping expansion of the Green 
function to higher orders [291 ]. Furthermore, a promis- 
ing way to obtain more insight into the superfluid phase 
would necessitate to combine the present technique with 
the effective action approach from Ref. [28| . 
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